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Abstract. It is shown using the Gross-Pitaevskii equation that resonance states 
of Bose-Einstein condensates with attractive interactions can be stabilized into 
true bound states. A semiclassical variational approximation and an independent 
quantum variational numerical method are used to calculate the energies (chemical 
potentials) and linewidths of resonances of the time-independent Gross-Pitaevskii 
equation; both methods produce similar results. Borders between the regimes 
of resonances, bound states, and, in two and three dimensions, collapse, are 
identified. 
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Resonances are metastable quantum states which have finite lifetimes; they play 
a key role in different types of scattering phenomena P^. Simon et al. [2] provided 
conditions for sharp transitions from a resonance to a true bound state as a parameter 
in the system's Hamiltonian is varied. These quantum phase transitions have been 
studied theoretically and experimentally in the context of the linear Schrodinger 
equation Our objective here is to search for similar transitions in nonlinear 

systems, such as Bosc-Einstcin condensates (BECs) 0]. In particular, we address 
the question of whether resonance states of a BEC with attractive binary atomic 
interactions can be stabilized into bound states by making the effective nonlinearity 
increasingly negative (through increase of the number of atoms, or tuning the s- 
wave scattering length a s near a Feshbach resonance to be more negative). As the 
existence and stability of eigenfunctions of nonlinear wave equations is also a central 
issue in other fields, and nonlinear Schrodinger type equations apply to many of 
them, including nonlinear optics, spin waves in magnetic films, Langmuir waves in 
hot plasmas, and gravity surface waves in fluids, the study of the effect of nonlinearity 
on resonances is of general physical interest. 

We develop two independent complementary approaches to this problem. The 
first utilizes a variational Gaussian ansatz for the quasi-stationary state and the WKB 
approximation to calculate the tunneling rate. The second utilizes a fully quantum 
variational method with an absorbing potential, which, in the case of the linear 
Schrodinger equation, may be shown formally [5] to be equivalent to the exterior 
scaling transformation [Jj]. We note parenthetically that complex scaling has been 
previously used to predict resonances for nonlinear Hartree-Fock models arising in 
electronic structure calculations of atoms and molecules 0. In both approaches, 
we seek to determine the energies and lincwidths of resonance states of the Gross- 
Pitaevskii equation (GPE), which is a nonlinear Schrodinger equation (NLS), derived 
in the mean-field approximation. In the context of the BEC, the "energy" may be 
interpreted as the chemical potential which is the real part of the eigenvalue associated 
with a resonance state, i.e., Re(/i). The "lincwidth" 7 may be interpreted as the 
coefficient in the rate equation, dN/dt = —-f(N)N, where N is the number of trapped 
atoms in the BEC and 7 (TV) = —2 (Im(/x)) /h is the decay rate. We emphasize that 
Im(/i) is a function of N: thus the tunneling problem associated with the decay of 
the resonance is a nonlinear one, and the lifetime of a resonance can be obtained by 
integrating dN/[—'y(N)N}. In the BEC, such nonlinear tunneling problems have been 
studied theoretically for bright solitons [Hj and experimentally in the context of the 
relaxation of condensate spin domains [§]. 

The two general methods are applied to the specific case of an external potential 
in the form of a harmonic well multiplied by a Gaussian envelope. This choice of the 
potential is motivated by studies of optically trapped BECs ^Oh where the waist of 
the trapping laser beam creates a broad Gaussian envelope. The central harmonic 
well may be created by a second, more narrowly focused laser beam, or by a higher 
Hermite-Gauss mode of the laser, as used in creating BEC waveguides [TI], or by a 
laser beam with embedded vorticity induced by passing the beam through a phase 
mask 12,. The lifetime given by tunneling of a quasi-bound resonant state can be 
reduced by other loss processes, such as three-body recombination and interactions 
with the background gas in the vacuum, as in the case for any trapped BEC. However, 
as these loss processes are the same for both resonant and bound states, tunneling of 
the mean field can cause a measurable reduction in lifetime. An external potential in 
the GPE is similar to the waveguiding profile created by a transverse modulation of the 
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refractive index that appears in the generalized NLS used for modeling propagation 
of spatial beams in optical media in effectively ID (planar waveguide) and 2D (bulk) 
geometries The potential that we introduce below corresponds to the waveguiding 
profile that can be easily engineered in optical media. 

The time-independent isotropic GPE with such an external potential can be 
written in terms of dimcnsionlcss variables scaled to the harmonic oscillator frequency 
uj and length ^ho = y/fi/rnuj as 

- + (D — l)tf'/r] /2 + V(r)tf + Uo\V 2 \V = , (1) 
V(r) = (V + r 2 /2) exp (-ar 2 ) , (2) 

(3d drr D - l \^>\ 2 = 1. (3) 
Jq 

Here ^(r) is the normalized wave function, p is the above-mentioned complex 
eigenvalue, /?i,2,3 = 2,27r,47r for the spatial dimension D = 1,2,3, and Uo oc a s N 
is the appropriately scaled nonlincarity coefficient in dimension D |14j . with a s the 
s-wave scattering length. The prime stands for d/dr, and a = (^ho/^Gauss) 2 is a 
parameter which characterizes the structure of the potential. For a broad Gaussian 
envelope, which pertains to the experimentally available optical trap described above, 
a < 1. 

Variational WKB approach — Approximate solutions to Eq. Q with potential 
can be obtained via variational methods. We adopt a Gaussian ansatz, ^ an (r) = 
A exp[— r 2 /(2p 2 )], where the width p and amplitude A are variational parameters. 
Substituting the ansatz into the normalization condition (J2J and the Lagrangian of 
the GPE (for the time being, Im(/x) is disregarded) and minimizing it with respect to 
p and A, one obtains A 2 = 2/[/3 D p D T(D/2)}, and 

Re(//) = {(~4 + D)/p 2 + {\ + a p 2 )- 2 - D ' 2 

x [Wo + (4 + D)p 2 - a(D + W Q a)p 4 } } /4 , (4) 

\U \ = T(D/2)2- 2+D / 2 f3p- 2+D (l + ap 2 )- 2 - D ' 2 

x {-2(1 + ap 2 ) 2+D ' 2 + /[2 - 4V a - a(D + 4V Q a)p 2 }} , (5) 

where T(D/2) is a Gamma function. The solution is stable in the framework of the 
time-dependent GPE if d~Rc(p) / d\Uo\ < 0, which is known as the Vakhitov-Kolokolov 
criterion |15j . 

One can now apply the WKB approximation to find lm(p) |16j . In the ID case, 
it is given by the standard expressions 

7 = exp(-2|/; i 2 drp(r)|)/T, (6) 
p = v/% - V cS (r)] , T = 4 | / ri dr/p(r)\ , (7) 

where V^ff(r) = \V(r) + Uo\*f?(r)\ 2 }, and the endpoints r = r\ and r — r2 arc found 
from p(r) = 0. In the 2D and 3D cases, it is necessary to transform the GPE into a 
ID form via 4* = 4>/\/r and \& = respectively (for the spherical wave). In the 
former case, this yields V e g(r) = V(r) — (l/8)r~ 2 + [/or _1 |0 2 |, whereas in the latter, 
V e g(r) = V(r) + U~or~ 2 \(f) 2 \. Note that lm(p) = —7/2 in these units; also, 7 must be 
multiplied by 2 in the ID case, to account for tunneling from both sides of the well. 
Results for the ID, 2D, and 3D cases are illustrated in the figures below. 

Exterior Complex Scaling of the GPE — In the linear Schrodingcr equation 
(Uq = 0), resonances are associated with complex eigenvalues of the Hamiltonian which 
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are obtained by imposing out going- wave boundary conditions on ^ . Such eigenvectors 
are sometimes called Siegert states. They are not in the Hilbert space, hence the 
Hamiltonian is a non-Hermitian operator over the class of functions including the 
Siegert states. In order to calculate resonances, it is convenient to carry out a 
similarity transformation which "brings back" the resonances to the Hilbert space. 
A well-known transformation of this kind is complex scaling, wherein the coordinates 
are rotated into the complex plane by a fixed angle 9 (El Q2| • When 9 is sufficiently 
large, the resonance eigenfunctions become square-integrable. Although the complex 
eigenvalues are 9- independent, provided tan(20) > 7/[2Re(/^)], the corresponding- 
resonance eigenfunctions do depend on 9: therefore the nonlinearity prevents direct 
application of complex scaling to the GPE. However, with the use of exterior |B] or 
smooth exterior scaling ^J|20|, the potential V(r) remains unsealed except at large r, 
and the small distortion of the effective nonlinear potential, [/ol^l 2 , may be neglected. 

It has been shown for the linear Schrodingcr equation that that smooth exterior 
scaling transformations are equivalent to adding a reflection-free complex absorbing 
potential (CAP) to the Hamiltonian [201, which is often approximated by a local 
negative imaginary potential. We chose Vcap(t) = —iX(r — ro) 4 , if r > ro, and 
VcapM = 0, if r < ro, where ro is a point near the edge of the grid, and A is a 
variational parameter (unrelated to the above variational approximation). This CAP 
is the most common one used in calculations of resonances in atoms, molecules, and 
nuclei. The value of A is found numerically from the variational condition d/i/dA = 0, 
where [i is generally complex, and stationary solutions with respect to small variation 
of A are obtained; note that stationary solutions in the A-variational space reguires 
optimization of another parameter, such as the box size in physical space, L (see 
the "cusp" condition in ^7]). The value of ro is chosen so that the wavefunction 
is exponentially small for r > ro; this ensures that the smooth exterior-scaling 
transformation docs not distort the effective potential V e g (r) (for a detailed discussion 
of this point, see 01201 and references therein). 

In our calculations, sinusoidal functions were used as a basis set in the ID and 3D 
cases, and Besscl functions in 2D. Starting from zero, the nonlinear strength parameter 
Uo was decreased in steps of SUq < 0, so that the overlap between the solutions for 
Uo and Uo + SUo was close to unity, to ensure that the evolution of a single resonance 
was followed. Thus, resonances could be identified by the CAP method, at each step 
of the iteration, in the same manner as in linear quantum mechanics. 

GPE Resonances and Bound States — Consider first the case of one dimension 
and no offset, i.e., D = 1 and V$ = 0. No bound states and a single resonance exist in 
this potential for the linear-Schrodinger (Uo = 0); the latter has Re(/i) = 0.41134 and 
Im(/i) = —0.0027894. Beginning with this known result, the Re(/z) and Im(/x) were 
calculated numerically as a function of the nonlinear parameter Uq by adiabatically 
decreasing the nonlinear parameter from zero. These numerical results, and those 
obtained via the above analytical methods, are displayed in Fig. ^ The transition 
from the resonance to a true bound state occurs at Uq 11 * = — 1.09. It is noteworthy 
that, despite the well known fact that any ID symmetric potential with minimum 
below the asymptotic value of the potential gives rise to at least one bound state 
in linear quantum mechanics, the nonlinear ID GPE creates a bound state only if 
-U > |C/ crit l (i.e., if C/ < C/ crit ). 

In the 2D case with Vb = 0, the numerical solution of Eqs. and JSJ shows that 
a resonance or bound state can exist for a < a c = 0.19245, where the tunneling rate 
diverges as a c — a — > +0. For a > a c , the potential walls arc too thin to retain the 
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condensate. The cases of a = 0.18 and a = 0.16 are shown in Fig. El The numerical 
study demonstrates that collapse occurs at the point Uq° 11 /(2it) = —0.86, while the 
variational approximation gives a result 13% larger (because the above ansatz is a poor 
approximation near the collapse, see Ref. [5] and references therein). The transition 
from the resonance to a bound state takes place at the point [/g nt /(2ir) = —0.66. 
These values vary slightly with a. 

In the 3D case, we again start by considering Vo = 0. Varying a, we found that 
resonances never turn into bound states as the nonlincarity coefficient Uo is made more 
negative. Instead, the condensate collapses as shown in Fig. El the variational 
approximation produces exactly the same result (details are not displayed here) . It is 
not surprising that collapse sets in after the stabilization of the bound state in the 2D 
case (Fig. |2Jl, while this does not happen in 3D: as is well known, the cubic nonlincarity 
gives rise to weak and strong collapse in 2D and 3D, respectively (the nonlinear term 
of the effective potential in 2D and 3D is Uor" 1 ]^] 2 and Uor~ 2 \cj)\ 2 , respectively). 

The above result for the 3D case may be understood as follows: the variational 
approximation with the Gaussian ansatz shows that the chemical potential at the 
collapse point in the isotropic 3D harmonic well without a Gaussian envelope is 
^coii _ o 224. Since this value is greater than zero, and the full potential has 
V(oo) — 0, it is not possible to transform a resonance into a bound state. However, 
lowering V(0) = Vo [see Eq. (J2J)] with respect to V(oo) makes such a transition possible. 
To illustrate this, the case of Vo = —0.8 and a = 0.02 is shown in Fig. 0] The 
corresponding potential (J2J) has a single resonance and no bound states in the linear 
limit (Uq = 0). The figure shows that, as \Uq\ increases, the width and energy of 
the resonance decrease, and indeed the resonance is stabilized into a bound state 
at the point U"'^ /{Air) = -0.552. As in the ID and 2D cases, the variational WKB 
approximation provides for reasonable accuracy, as compared to the numerical results. 

Clearly, our two independent methods may also be used to address the question of 
whether BECs with repulsive nonlincarity have a critical point beyond which bound 
states trapped in a potential well destabilize into resonances as the nonlincarity is 
increased. 

In conclusion, we have demonstrated within a mean field approximation that 
resonances in a BEC can be transformed into true bound states, as the strength of the 
attractive nonlinearity increases beyond a critical value |[/g rit |. Borders between three 
different dynamical regimes, viz., resonances, bound states, and collapse, have been 
delineated. In ID, resonances can always be stabilized into bound states, as collapse 
is not possible. In 2D, we find that the transition of the resonance into a bound state 
can be readily achieved before collapse occurs, as the collapse is weak in this case. 
However, in 3D, where collapse is strong, the critical value for the transition, |t/o nt |j 
can be made smaller than the collapse threshold, | C/q o11 | , only if V(0) < V(oo), i.e., the 
floor of the potential well is lower than at infinity. In optically trapped BECs, this can 
be achieved by means of depression of the effective potential well with a red-detuned 
laser. 
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Figure 1. One dimension. Shown is the Re(/i) and Im(/i) versus | fTo I for 
the potential in the form of a harmonic well times the Gaussian envelope [see 
Eq. with a = (^ho/^Gauss) 2 = 0-2. Solid curves: results of the numerical 

method utilizing a complex scaling method. Dashed curves: the variational WKB 
approximation. The critical point for the conversion of the resonance into a bound 
state is \U^ it \ = 1.09. 
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Figure 2. Two dimensions. Same as in Fig. Q for two different values of 
the potential-shape parameter, a = 0.16 and a = 0.18 (the upper and lower 
curves, respectively; the analytical curves for Re(/i) at both values of a completely 
overlap). Regions of the resonance, bound state, and collapse are indicated. 
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Figure 3. Three dimensions. The width of the resonance states versus energy for 
three different wells with V(0) = V(oo) (i.e., V =0): a = (^ ho /f Gauss) 2 = 0.02, 
0.03, and 0.04. For the definition of k see the label attached to the vertical axis. 
In each case, the collapse point is reached before the resonance can be stabilized 
into a bound state. The variational WKB approximation produces similar results 
(not shown here). 




Figure 4. Three dimensions. Same as in Fig.lTlfor the potential J2J with a = 0.02 
and Vo = -0.8 (so that V(0) = V(oo) - 0.8). The offset Vb < allows for the 
stabilization of the resonance into a bound state, unlike the case shown in Fig. 151 



